# simplified version of StrategyUnitTheme::su_theme_pal()
#   see: https://github.com/The-Strategy-Unit/StrategyUnitTheme
su_theme_pal <- function(palette = "main", ...) {
  pal <- c(orange = "#f9bf07",
           red    = "#ec6555",
           blue   = "#5881c1")
  grDevices::colorRampPalette(pal, ...)
}

If you haven’t already got these packages installed you will need the following to complete this workshop:

install.packages(c("tidyverse", "simmer", "simmer.plot", "gridExtra", "TruncatedNormal"))

1 Building the model - part 1

1.1 Libraries needed

suppressPackageStartupMessages({
  library(tidyverse)
  library(simmer)
  library(simmer.plot)
  library(TruncatedNormal)
  library(writexl)
  library(gridExtra)
})

1.2 Create distribution functions

For the patients interarrival rate, the time taken for asking the screening questions and the time for the receptionists to make the decision on what appointment to book.

You can use helper functions

# generic truncated normal distributions looks like:
rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)

## [1] 0.6213012

dist_helper <- function(rfn, ...) {
  function() rfn(1, ...)
}
dist_patient_interarrival <- dist_helper(rexp, 50 / 60)
dist_screening            <- dist_helper(rtnorm, 4, 1,    0, Inf)
dist_decision             <- dist_helper(rtnorm, 1, 0.25, 0, Inf)
dist_booking              <- dist_helper(rtnorm, 1, 0.25, 0, Inf)

Alternatively:

dist_patient_interarrival <- function() rexp(1, 50 / 60)
dist_screening            <- function() rtnorm(n = 1, mu = 4, sd = 1,    lb = 0, ub = Inf)
dist_decision             <- function() rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)
dist_booking              <- function() rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)
dist_contact_mode         <- function() sample(1:2, 1, FALSE, c(0.15, 0.85))
dist_booking_type         <- function() sample(1:3, 1, FALSE, c(0.3, 0.6, 0.1))

1.3 Creating the trajectories for the branches

Think about the flow diagram of the process to create the trajectories for the model, let’s break it down into manageable sections.

1.3.1 ECR branch

Create a trajectory for this branch called branch_ecr. Within the trajectory use set_attribute via_ecr to label the patient with a value of 1, so we will know what route they have come to analyse later. Then seize a receptionist ready for the next stage of the process.

branch_ecr <- trajectory("branch ecr") %>%
  set_attribute("via_ecr", 1) %>%
  seize("receptionist")
branch_ecr

## trajectory: branch ecr, 2 activities
## { Activity: SetAttribute | keys: [via_ecr], values: [1], global: 0, mod: N, init: 0 }
## { Activity: Seize        | resource: receptionist, amount: 1 }

1.3.2 Phone branch

create a trajectory for this branch called branch_phone. Within the trajectory use set_attribute via_ecr to label the patient with a value of 0 (as opposed to 1 which are ECR route patients).

The patients renege after 5 minutes so log this and set the attribute "reneged" to 1. Then seize a receptionist ready for the next stage. Once the patient has a receptionist abort the renege and add timeout of distribution screening

renege_hungup <- trajectory("hung up") %>%
  set_attribute("reneged", 1) %>%
  log_("lost my patience. Hanging up...")

branch_phone <- trajectory("branch phone") %>%
  set_attribute("via_ecr", 0) %>%
  renege_in(t = 5, out = renege_hungup) %>%
  seize("receptionist") %>%
  renege_abort() %>% # now have a receptionist the patient won't renege
  timeout(dist_screening)

1.3.3 Booking option branches

Receptionist is deciding if a patient needs phone appointment, needs phone appointment or does not need any appointment.

branch_booking_none <- trajectory("no appointment needed") %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 0)
branch_booking_none

## trajectory: no appointment needed, 2 activities
## { Activity: Release      | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [0], global: 0, mod: N, init: 0 }

branch_booking_phone <- trajectory("booking a phone") %>%
  timeout(dist_booking) %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 1)
branch_booking_phone

## trajectory: booking a phone, 3 activities
## { Activity: Timeout      | delay: function() }
## { Activity: Release      | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [1], global: 0, mod: N, init: 0 }

branch_booking_f2f <- trajectory("booking a f2f") %>%
  timeout(dist_booking) %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 2)
branch_booking_f2f

## trajectory: booking a f2f, 3 activities
## { Activity: Timeout      | delay: function() }
## { Activity: Release      | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [2], global: 0, mod: N, init: 0 }

1.4 Patient Trajectory

patient_flow <- trajectory("patient flow") %>%
  branch(dist_contact_mode, TRUE,
         branch_ecr,
         branch_phone) %>%
  timeout(dist_decision) %>%
  branch(dist_booking_type, TRUE,
         branch_booking_none,
         branch_booking_phone,
         branch_booking_f2f)
plot(patient_flow, fill = su_theme_pal("main"))

# you can change the verbose argument to TRUE to get more information
plot(patient_flow, fill = su_theme_pal("main"), verbose = TRUE)

1.5 Add resources and generator to the environment

You can use to() function in the generator, so that it completes the arrivals. Use to() with caution as the stop time is that of the generator and not the simulation. The simulation will run until all arrivals are finished, so we are assuming that staff stay overtime

env_part1 <- simmer("part 1") %>%
  add_resource("receptionist", 2) %>%
  add_generator("patient", patient_flow, to(8 * 60, dist_patient_interarrival), mon = 2)

1.6 Run the model

1.7 Glimpse results to check run

env_part1 %>%
  get_mon_arrivals() %>%
  arrange(start_time) %>%
  head(10)

name	start_time	end_time	activity_time	finished	replication
patient0	0.9062182	8.666809	7.760591	TRUE	1
patient1	1.0810663	2.880660	1.799594	TRUE	1
patient2	1.6043486	7.786732	4.906072	TRUE	1
patient3	2.0532035	7.053204	0.000000	TRUE	1
patient4	2.3693078	7.369308	0.000000	TRUE	1
patient5	2.7333249	7.733325	0.000000	TRUE	1
patient6	3.6477607	13.937513	6.150780	TRUE	1
patient7	3.9308550	13.606157	4.939348	TRUE	1
patient8	4.9175231	9.917523	0.000000	TRUE	1
patient9	7.7549414	12.754941	0.000000	TRUE	1

env_part1 %>%
  get_mon_arrivals() %>%
  arrange(start_time) %>%
  tail(10)

	name	start_time	end_time	activity_time	finished	replication
389	patient388	472.4713	477.9166	5.445350	TRUE	1
390	patient389	472.8399	483.2193	6.273473	TRUE	1
391	patient390	473.1776	481.5829	3.666322	TRUE	1
392	patient391	474.4332	479.4332	0.000000	TRUE	1
393	patient392	475.5632	480.5632	0.000000	TRUE	1
394	patient393	476.6170	487.7046	6.121711	TRUE	1
395	patient394	476.7175	481.7175	0.000000	TRUE	1
396	patient395	477.5498	482.5498	0.000000	TRUE	1
397	patient396	478.8235	490.1992	6.979876	TRUE	1
398	patient397	479.2955	484.2955	0.000000	TRUE	1

1.8 Questions

If you are new to R and struggling with data manipulation, you can export the results as xls or csv. Or email us and we can help you!

write_xlsx(get_mon_arrivals(env_part1), "arrivals_part1.xlsx")
write_xlsx(get_mon_resources(env_part1), "resources_part1.xlsx")
write_xlsx(get_mon_attributes(env_part1), "attributes_part1.xlsx")

1.8.1 Count reneged

env_part1 %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  summarise(countreneged = n())

countreneged
179

n = 148 with 2 receptionists

1.8.2 Count F2F

env_part1 %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  summarise(countF2F = n())

countF2F
30

n = 27 with 2 receptionists

1.8.3 Call waits

Get attributes of finished patients from get_mon_attributes(), join with get_mon_arrivals(). Because of varying number of decimal places we get some rounding errors resulting in near zero but negative results. That’s why when we are doing data wrangling, wee use pmax (), so it takes 0 or the wait time which ever is greater. We can then calculate summary statistics.

part1_attributes <- env_part1 %>%
  get_mon_attributes() %>%
  dplyr::select(-time) %>%
  pivot_wider(names_from = "key", values_from = "value", values_fill = 0)

part1_callwait_results <- env_part1 %>%
  get_mon_arrivals(ongoing = TRUE) %>% # not needed as there aren't any due to the generator and run()
  inner_join(part1_attributes, by = c("name", "replication")) %>%
  filter(via_ecr == 0) %>%
  mutate(wait_time_call = case_when(
    # the caller reneged
    reneged == 1 ~ 5,
    # in queue still - these could also still be in progress
    # (would need to know time resource seized to calculate wait time),
    # except we used to() in generator and run() so we have no unfinished
    # !finished ~ (10 * 60) - start_time,
    # everything else, (assume has been answered)
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

part1_callwait_summary <- part1_callwait_results %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished))

part1_callwait_summary

median_wait	mean_wait	var_wait	min_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countunfinished	countfinished
5	4.223388	1.733617	0	5	3.906749	5	0	329

with 2 receptionists

#     median_wait   mean_wait      var_wait    min_wait    max_wait      Qtile25th_wait    Qtile75th_wait
#1    4.839153        3.811972      2.657451         0         5           3.010095              5

#     countunfinished countfinished
#1               0           334

2 Scenario additional model features

Glimpse the end tail and also check whether ECR backlog is cleared at the end of the days. Overall backlog - let’s look at 1 replication

attributes <- env_part1 %>%
  get_mon_attributes() %>%
  filter(replication == 1)

test <- env_part1 %>%
  get_mon_arrivals() %>%
  filter(finished == FALSE, replication == 1) %>%
  inner_join(attributes, by = "name") %>%
  arrange(start_time)

test

name	start_time	end_time	activity_time	finished	replication.x	time	key	value	replication.y

env_part1 %>%
  get_mon_attributes() %>%
  filter(replication == 1, key == "via_ecr", value == 1) %>%
  count()

n
69

You can see that everyone finished and we had total of 66 ECR.

3 Extra model features - part 2

3.1 New environment

Create a new environment

env_part2 <- simmer("part2")
env_part2

## simmer environment: part2 | now: 0 | next: 
## { Monitor: in memory }

3.2 Add time-dependent interarrival times

To include time dependent interarrivals create a different distribution, one for each time period (3 in total). You may also find using these helper functions to switch between minutes, hours and days useful:

t_minute <- function(t) t
t_hour   <- function(t) t_minute(t) * 60
t_day    <- function(t) t_hour(t) * 24

dist_patient_interarrival1 <- function() rexp(1, 120 / t_hour(2))
dist_patient_interarrival2 <- function() rexp(1, 240 / t_hour(6))
dist_patient_interarrival3 <- function() rexp(1,  40 / t_hour(2))

3.3 Phone branch with prioritisation

We are also going prioritise the patients contacting the surgery via the phone.

In the ECR branch, no changes are needed. The default priority for all arrivals is set up to 0.

In the phone branch, we want to set a higher priority using set_prioritization(), that receptionists prioritise answering calls over reviewing electronic requests but without interrupting current tasks (prioirity = 1, preemptive = 1, restart = FALSE). If you didn’t use now(env_part1) to set attributes you can use join(branch_phone) instead to link after the prioritisation back to your original phone branch.

branch_phone2 <- trajectory("branch phone prioritised") %>%
  set_prioritization(c(1, 1, FALSE)) %>%
  join(branch_phone)

plot(branch_phone2, fill = su_theme_pal("main"))

from the help page for add_generator():

priority: the priority of each arrival (a higher integer equals higher priority; defaults to the minimum priority, which is 0.
preemptible: if a seize occurs in a preemptive resource, this parameter establishes the minimum incoming priority that can preempt these arrivals (an arrival with a priority greater than preemptible gains the resource). In any case, preemptible must be equal or greater than priority, and thus only higher priority arrivals can trigger preemption.
restart: whether the activity must be restarted after being preempted.

3.4 Patient Trajectory

same as Part 1 but with the new prioritised phone branch

patient_flow2 <- trajectory("patient flow2") %>%
  branch(dist_contact_mode, TRUE,
         branch_ecr,
         branch_phone2) %>% #<-- changed the phone branch
  timeout(dist_decision) %>%
  branch(dist_booking_type, TRUE,
         branch_booking_none,
         branch_booking_phone,
         branch_booking_f2f)
plot(patient_flow2, fill = su_theme_pal())

3.5 Schedule for staffing

Use the schedule(timetable, values, period) function and set it to repeat the schedule every 24 hour period. Work out from the numbers on each shift what times capacity changes and how many are working at that point:

e.g. 0-8am = 0, 8-9am = 1, 9-9.30 am etc. becomes:

sch_receptionists <- schedule(
  timetable = t_hour(c(0.0, 8.0, 9.0, 9.5, 12.0, 15.5, 18.0)),
  values    =        c(0,   1,   2,   4,    3,    1,    0),
  period    = t_hour(24)
)

3.6 Add resources and generator to the environment

Add the receptionists as a resource using the schedule you created for the capacity. Either add three different generators to the model environment setting the times they operate with from_to or use functional programming to control this. You can look back at the training slides and examples to help you.

env_part2 %>%
  add_resource("receptionist", capacity = sch_receptionists) %>%
  add_generator("patient[8am-10am]", patient_flow2,
                from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
                mon = 2) %>%
  add_generator("patient[10am-4pm]", patient_flow2,
                from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
                mon = 2) %>%
  add_generator("patient[4pm-6pm]", patient_flow2,
                from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
                mon = 2)

## simmer environment: part2 | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: receptionist | monitored: TRUE | server status: 0(0) | queue status: 0(Inf) }
## { Source: patient[8am-10am] | monitored: 2 | n_generated: 0 }
## { Source: patient[10am-4pm] | monitored: 2 | n_generated: 0 }
## { Source: patient[4pm-6pm] | monitored: 2 | n_generated: 0 }

3.7 Run model

Set the random seed to 1 to give us the same results each time. Run your model for 5 days (52460) or use one of the time helper functions)

3.8 Questions

Repeat your analysis of call waits, reneging and F2F appointments as you did for part 1 using the environment for part 2. Plot utilisation of the resources using the plot function with get_mon_arrivals() and the metric “waiting_time”.

3.8.1 Count reneged

env_part2 %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  summarise(countreneged = n())

countreneged
582

n = 582 with 2 receptionists

3.8.2 Count F2F

env_part2 %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  summarise(countF2F = n())

countF2F
160

n = 160 with 2 receptionists

3.8.3 Call waits

part2_attributes <- env_part2 %>%
  get_mon_attributes() %>%
  dplyr::select(-time) %>%
  pivot_wider(names_from = "key", values_from = "value", values_fill = 0)

part2_callwait_results <- env_part2 %>%
  get_mon_arrivals(ongoing = TRUE) %>%
  full_join(part2_attributes, by = c("name", "replication")) %>% # or can use inner_join
  filter(via_ecr == 0) %>%
  mutate(wait_time_call = case_when(
    reneged == 1 ~ 5,
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

part2_callwait_summary <- part2_callwait_results %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished))

part2_callwait_summary

median_wait	mean_wait	var_wait	min_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countunfinished	countfinished
3.825883	3.183036	3.479728	0	5	1.464796	5	0	1726

#   median_wait     mean_wait     var_wait    min_wait    max_wait    Qtile25th_wait    Qtile75th_wait
#1    3.825883      3.183036      3.479728          0          5          1.464796              5
#   countunfinished countfinished
#1               0          1726

3.9 What else would you do to improve this model?

Some ideas:

prevent staff staying to finish all just the one they are on,
multiple runs,
staff breaks, availability, are they also doing other tasks?

3.10 What other metrics would you want to measure / the practice find useful?

Below are some ideas on additional analysis you might want to do

Plot waiting time over time

plot(get_mon_arrivals(env_part2), metric = "waiting_time")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

boxplot to show variation in metrics e.g. call_wait

g1 <- ggplot() +
  geom_boxplot(data = part2_callwait_results,
               aes(x = "rep=1", y = wait_time_call),
               outlier.shape = 19,
               width = 0.3,
               fill = "white",
               colour = "orange") +
  geom_point(data = part2_callwait_summary,
             aes(x = "rep=1", y = mean_wait),
             shape = 19,
             size = 2,
             colour = "red") +
  labs(x = "Part2",
       y = "call wait (mins)",
       title = "Boxplot of call wait from single replication") +
  theme(axis.text = element_text(size = 11))

g1

Glimpse the end tail and also check whether ECR backlog is cleared at the end of the days. Overall backlog:

attributes <- env_part2 %>%
  get_mon_attributes() %>%
  filter(replication == 1)

test <- env_part2 %>%
  get_mon_arrivals() %>%
  filter(finished == FALSE, replication == 1) %>%
  inner_join(attributes, by = "name") %>%
  arrange(start_time)

test %>%
  filter(key == "via_ecr", value != 1)

name	start_time	end_time	activity_time	finished	replication.x	time	key	value	replication.y

env_part2 %>%
  get_mon_attributes() %>%
  filter(replication == 1, key == "via_ecr", value == 1) %>%
  count()

n
357

test %>%
  head(10) %>%
  as_tibble()

name	start_time	end_time	activity_time	finished	replication.x	time	key	value	replication.y
patient[4pm-6pm]191	6774.405	NA	NA	FALSE	1	6774.405	via_ecr	1	1
patient[4pm-6pm]192	6775.585	NA	NA	FALSE	1	6775.585	via_ecr	1	1
patient[4pm-6pm]201	6790.286	NA	NA	FALSE	1	6790.286	via_ecr	1	1

All unfinished are ECR, because calls were prioritised ahead of them. 357 total ECR.

# A tibble: 3 x 10 - 3 are unfinished in a week out of 357
# name                start_time end_time activity_time finished replication.x  time key     value replication.y
#<chr>                    <dbl>    <dbl>         <dbl> <lgl>            <int> <dbl> <chr>   <dbl>         <int>
#  patient[4pm-6pm]191      6774.       NA            NA FALSE                1 6774. via_ecr     1             1
#  patient[4pm-6pm]192      6776.       NA            NA FALSE                1 6776. via_ecr     1             1
#  patient[4pm-6pm]201      6790.       NA            NA FALSE                1 6790. via_ecr     1             1

Look at resources (reception staff in this example)

resources2 <- get_mon_resources(env_part2)
plot(resources2, metric = "usage", "receptionist", items = "server", steps = TRUE)

plot(resources2, metric = "utilization", "receptionist")

# show the top 10 results
head(resources2, n = 10)

resource	time	server	queue	capacity	queue_size	system	limit	replication
receptionist	480.0000	0	0	1	Inf	0	Inf	1
receptionist	480.0000	1	0	1	Inf	1	Inf	1
receptionist	481.1816	1	1	1	Inf	2	Inf	1
receptionist	481.6397	1	2	1	Inf	3	Inf	1
receptionist	484.5346	1	3	1	Inf	4	Inf	1
receptionist	484.9699	1	4	1	Inf	5	Inf	1
receptionist	485.1100	1	3	1	Inf	4	Inf	1
receptionist	486.1832	1	4	1	Inf	5	Inf	1
receptionist	486.6397	1	3	1	Inf	4	Inf	1
receptionist	487.5739	1	4	1	Inf	5	Inf	1

# list of hours
HrList <- data.frame(
  resource = "receptionist",
  time = seq(from = 0, to = t_day(5), by = 60),
  server = NA,
  queue = NA,
  capacity = NA,
  queue_size = Inf,
  system = NA,
  limit = Inf,
  replication = 1
)

# which aren't in the data
AddSnaps <- filter(HrList, !(time %in% resources2$time))

resources2_hourly <- bind_rows(resources2, AddSnaps) %>%
  arrange(time) %>%
  mutate(server = ifelse(is.na(server), lag(server), server),
         queue = ifelse(is.na(queue), lag(queue), queue),
         capacity = ifelse(is.na(capacity), lag(capacity), capacity),
         system = ifelse(is.na(system), lag(system), system))

resourceutilisation_byhour <- resources2_hourly %>%
  arrange(time) %>%
  mutate(slot_hour = floor(time / 60)) %>%
  mutate(slot_day = ceiling(slot_hour / 24)) %>%
  mutate(slot = slot_hour - ((slot_day - 1) * 24)) %>%
  group_by(slot) %>%
  mutate(dt = time - lag(time),
         capacity = ifelse(capacity < server, server, capacity),
         in_use = dt * lag(server / capacity)) %>%
  summarise(utilised_percent = sum(in_use, na.rm = TRUE) / sum(dt, na.rm = TRUE) * 100)

resourceutilisation_byhour

slot	utilised_percent
1	0.0000000
2	0.0000000
3	0.0000000
4	0.0000000
5	0.0000000
6	0.0000000
7	0.0000000
8	100.0000000
9	100.0000000
10	99.9259080
11	63.7850307
12	99.6780740
13	91.7546945
14	75.7575154
15	99.8493045
16	99.8748793
17	100.0000000
18	0.1985734
19	0.0000000
20	0.0000000
21	0.0000000
22	0.0000000
23	0.0000000
24	0.0000000

3.11 Alternative way to calculate waiting time - option B

Option B. If you had used now(env_part1) to set attributes then you will need to copy your code and change the environment for all parts of the trajectory in which it was used. You then need to create the overall trajectory again using the new version of the phone branch (and other parts if you have used now()).

Firstly, create a new environment.

env_part2B <- simmer("part 2B")
env_part2B

## simmer environment: part 2B | now: 0 | next: 
## { Monitor: in memory }

Change references to part env_part2B in all trajectories

renege_hungup2B <- trajectory("hung up") %>%
  set_attribute("reneged", 1) %>%
  set_attribute("reneged_at", function() now(env_part2B)) %>% # changed to env_part2B
  log_("lost my patience. Hanging up...")

branch_phone2B <- trajectory("branch phone") %>%
  set_prioritization(c(1, 1, FALSE)) %>% # set to increase priority from default of 0 to 1, no preempt and no restart
  set_attribute("via_ecr", 0) %>%
  set_attribute("arrival", function() now(env_part2B)) %>% # changed to env_part2B
  renege_in(t = 5, out = renege_hungup) %>%
  seize("receptionist") %>%
  set_attribute("call_taken", function() now(env_part2B)) %>% # changed to env_part2B
  renege_abort() %>%
  timeout(dist_screening)

branch_booking_none2B <- trajectory("no appointment needed") %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 0) %>%
  set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B

branch_booking_phone2B <- trajectory("booking a f2f") %>%
  timeout(dist_booking) %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 1) %>%
  set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B

branch_booking_f2f2B <- trajectory("booking a f2f") %>%
  timeout(dist_booking) %>%
  release("receptionist") %>%
  set_attribute("booked_appt", 2) %>%
  set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B

patient_flow2B <- trajectory("patient flow2B") %>%
  branch(dist_contact_mode, TRUE,
         branch_ecr,
         branch_phone2B) %>%
  timeout(dist_decision) %>%
  branch(dist_booking_type, TRUE,
         branch_booking_none2B,
         branch_booking_phone2B,
         branch_booking_f2f2B)

Populate environment using the trajectory for patient flow 2B

env_part2B %>%
  add_resource("receptionist", capacity = sch_receptionists) %>%
  add_generator("patient[8am-10am]", patient_flow2B,
                from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
                mon = 2) %>%
  add_generator("patient[10am-4pm]", patient_flow2B,
                from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
                mon = 2) %>%
  add_generator("patient[4pm-6pm]", patient_flow2B,
                from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
                mon = 2)

## simmer environment: part 2B | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: receptionist | monitored: TRUE | server status: 0(0) | queue status: 0(Inf) }
## { Source: patient[8am-10am] | monitored: 2 | n_generated: 0 }
## { Source: patient[10am-4pm] | monitored: 2 | n_generated: 0 }
## { Source: patient[4pm-6pm] | monitored: 2 | n_generated: 0 }

Run your model for 5 days

3.11.1 Part 2 Questions for Option B

What is the average wait for calls to be answered, number abandoning their calls (hanging up) and the number booked for a face to face appointment in this simulated week?

env_part2B %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  summarise(countreneged = n())

countreneged
582

n=582

How many people get a F2F booking?

env_part2B %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  summarise(countF2F = n())

countF2F
160

n=160

What are the call waits?

part2B_callwait_results <- env_part2B %>%
  get_mon_attributes() %>%
  pivot_wider(-c(time, replication), names_from = key, values_from = value) %>%
  filter(via_ecr == 0, arrival >= 0) %>%
  mutate(wait_time_call = case_when(
    !is.na(reneged) ~ 5,
    is.na(call_taken) ~ (10 * 60) - arrival, # in queue/progress still, model run for 10 hours
    TRUE ~ call_taken - arrival # has been answered
  ))

part2B_callwait_summary <- part2B_callwait_results %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(is.na(call_end) & is.na(reneged)),
            countfinished = sum(!is.na(call_end) | !is.na(reneged)))

part2B_callwait_summary

median_wait	mean_wait	var_wait	min_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countunfinished	countfinished
3.825883	3.183036	3.479728	0	5	1.464796	5	0	1726

#   # A tibble: 1 x 9
#       median_wait     mean_wait     var_wait  min_wait    max_wait Qtile25th_wait Qtile75th_wait
#1         3.83          3.18        3.48        0              5           1.46              5
# countunfinished countfinished
# 0                 1726

4 Multiple runs and Scenarios - part 3

Now we are going to run our model 50 times and analyse different scenarios.

4.1 Replication of baseline

Use lapply to loop through and replicate the environment you created in Part 2 100 times

4.2 get baseline results and plots

Update some of the code you used to look at results in pasrt 1 and 2. This time group_by(replication) before you summarise. E.g count of the reneged, then ungroup(). This will give the result for each replication. It is these that are then summarised again to give your final results for i.e. the mean of the mean number reneged in each replication.

4.2.1 Count reneged

envs_baseline %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  group_by(replication) %>%
  summarise(countreneged = n()) %>%
  ungroup() %>%
  summarise(countreneged_mean = mean(countreneged))

countreneged_mean
579.86

n=580

4.2.2 Count F2F

envs_baseline %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  group_by(replication) %>%
  summarise(countF2F = n()) %>%
  ungroup() %>%
  summarise(countF2F_mean = mean(countF2F))

countF2F_mean
143.04

n=143

Before we look at call waits let’s understand a little more about unfinished. using the staff schedule means that the final staff member finishes at 6pm on the 5th day so at (42460)+(18*60)= 6840, there are 4 in the system, 1 patient in progress and 3 waiting the staff member finishes that one at 6845 and then leaves with 3 remaining unfinished a similar thing occurs at the end of day 3:

get_mon_resources(envs_baseline) %>%
  filter(replication == 1) %>%
  tail(n = 10) # show last 10 rows

	resource	time	server	queue	capacity	queue_size	system	limit	replication
4184	receptionist	6824.155	1	4	1	Inf	5	Inf	1
4185	receptionist	6827.207	1	5	1	Inf	6	Inf	1
4186	receptionist	6829.155	1	4	1	Inf	5	Inf	1
4187	receptionist	6829.959	1	3	1	Inf	4	Inf	1
4188	receptionist	6832.803	1	4	1	Inf	5	Inf	1
4189	receptionist	6835.244	1	3	1	Inf	4	Inf	1
4190	receptionist	6835.407	1	4	1	Inf	5	Inf	1
4191	receptionist	6839.038	1	3	1	Inf	4	Inf	1
4192	receptionist	6840.000	1	3	0	Inf	4	Inf	1
4193	receptionist	6845.235	0	3	0	Inf	3	Inf	1

get_mon_resources(envs_baseline) %>%
  filter(replication == 1, time < 6240) %>%
  tail(n = 10)

	resource	time	server	queue	capacity	queue_size	system	limit	replication
3326	receptionist	5387.154	1	4	1	Inf	5	Inf	1
3327	receptionist	5387.827	1	3	1	Inf	4	Inf	1
3328	receptionist	5391.761	1	4	1	Inf	5	Inf	1
3329	receptionist	5392.154	1	3	1	Inf	4	Inf	1
3330	receptionist	5392.374	1	4	1	Inf	5	Inf	1
3331	receptionist	5392.431	1	5	1	Inf	6	Inf	1
3332	receptionist	5393.999	1	4	1	Inf	5	Inf	1
3333	receptionist	5397.431	1	3	1	Inf	4	Inf	1
3334	receptionist	5400.000	1	3	0	Inf	4	Inf	1
3335	receptionist	5401.615	0	3	0	Inf	3	Inf	1

envs_baseline %>%
  get_mon_arrivals(ongoing = TRUE) %>% # 101972 rows
  filter(replication == 1, is.na(end_time))

name	start_time	end_time	activity_time	finished	replication
patient[4pm-6pm]212	-1.000	NA	NA	FALSE	1
patient[4pm-6pm]201	6790.286	NA	NA	FALSE	1
patient[4pm-6pm]191	6774.405	NA	NA	FALSE	1
patient[4pm-6pm]192	6775.585	NA	NA	FALSE	1
patient[10am-4pm]1218	-1.000	NA	NA	FALSE	1
patient[8am-10am]653	-1.000	NA	NA	FALSE	1

4.2.3 Call waits

baseline_attributes <- envs_baseline %>%
  get_mon_attributes() %>% #202732 obs
  dplyr::select(-time) %>%
  pivot_wider(names_from = "key", values_from = "value", values_fill = 0)
#101468 obs

baseline_callwait_results <- envs_baseline %>%
  get_mon_arrivals(ongoing = TRUE) %>%
  distinct() %>% # goes down to 101618, removes duplicates of unfinished
  # still more in arrivals than attributes as attributes only has finished
  full_join(baseline_attributes, by = c("name", "replication")) %>%
  filter(via_ecr == 0) %>% # phone only also removes those yet to start e.g. start_time -1
  mutate(wait_time_call = case_when(
    # the caller reneged
    reneged == 1 ~ 5,
    # unfinished still in the system,
    !finished ~ t_day(5) - start_time,
    # everything else, (assume has been answered)
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

## calculate summary statistics
baseline_callwait_replication_summary <- baseline_callwait_results %>%
  group_by(replication) %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished)) %>%
  mutate(scenario = "baseline")

# show top 10 rows
head(baseline_callwait_replication_summary, n = 10)

replication	median_wait	mean_wait	var_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countfinished	scenario
1	3.825883	3.183036	3.479728	5	1.464796	5	1726	baseline
2	3.744117	3.221769	3.258414	5	1.555424	5	1770	baseline
3	3.775021	3.116034	3.682359	5	1.206570	5	1772	baseline
4	3.952666	3.235613	3.446565	5	1.535425	5	1722	baseline
5	3.810451	3.201863	3.433046	5	1.446567	5	1736	baseline
6	3.550199	3.053321	3.688564	5	1.141632	5	1699	baseline
7	3.525437	3.005393	3.874786	5	1.005038	5	1652	baseline
8	3.808467	3.141373	3.616938	5	1.237069	5	1733	baseline
9	3.892609	3.168985	3.635974	5	1.225592	5	1750	baseline
10	3.623343	3.068928	3.714392	5	1.107349	5	1662	baseline

baseline_callwait_summary <- baseline_callwait_replication_summary %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished)) %>%
  mutate(scenario = "baseline") # to identify it later

baseline_callwait_summary

median_wait	mean_wait	var_wait	min_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countunfinished	countfinished	scenario
3.765179	3.147497	3.580269	0	5	1.334405	5	0	1722.8	baseline

5 Scenario

If a neighbouring practice closed leading to an increase in demand of 45%. Whilst the practice were able to better advertise the ECR route to patients such that 40% of all the patients (old and new) were able to use this route. What impact would this have on the call wait and reneging of patients?

5.1 create patient flow 3 for this scenario with 40% ECR

dist_contact_mode_S <- function() sample(1:2, 1, FALSE, c(0.4, 0.6))
patient_flow3 <- trajectory("patient flow3") %>%
  branch(dist_contact_mode_S, TRUE, # <- changed to the scenario distribution
         branch_ecr,
         branch_phone2) %>%
  timeout(dist_decision) %>%
  branch(dist_booking_type, TRUE,
         branch_booking_none,
         branch_booking_phone,
         branch_booking_f2f)

5.2 scenario interarrivals

dist_patient_interarrival1_S <- function() rexp(1, (120 * 1.45) / t_hour(2))
dist_patient_interarrival2_S <- function() rexp(1, (240 * 1.45) / t_hour(6))
dist_patient_interarrival3_S <- function() rexp(1,  (40 * 1.45) / t_hour(2))

5.3 lapply as before but using patient_flow3 and increasing frquency of arrivals

5.4 Scenario results

Repurpose some of the code you used to look at results in parts 1 and 2 of the workshop. Remember to store the results with a different name to your baseline results.

5.4.1 Count reneged

envs_scenario %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  group_by(replication) %>%
  summarise(countreneged = n()) %>%
  ungroup() %>%
  summarise(countreneged_mean = mean(countreneged))

countreneged_mean
603.16

n=752

5.4.2 Count F2F

envs_scenario %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  group_by(replication) %>%
  summarise(countF2F = n()) %>%
  ungroup() %>%
  summarise(countF2F_mean = mean(countF2F))

countF2F_mean
171.56

n=143

5.4.3 Call waits

scenario_attributes <- envs_scenario %>%
  get_mon_attributes() %>%
  dplyr::select(-time) %>%
  pivot_wider(names_from = "key", values_from = "value", values_fill = 0)

scenario_callwait_results <- envs_scenario %>%
  get_mon_arrivals(ongoing = TRUE) %>%
  distinct() %>%
  full_join(scenario_attributes, by = c("name", "replication")) %>%
  filter(via_ecr == 0) %>%
  mutate(wait_time_call = case_when(
    # the caller reneged
    reneged == 1 ~ 5,
    # unfinished still in the system,
    !finished ~ t_day(5) - start_time,
    # everything else, (assume has been answered)
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

## calculate summary statistics
scenario_callwait_replication_summary <- scenario_callwait_results %>%
  group_by(replication) %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished)) %>%
  mutate(scenario = "scenario")

scenario_callwait_summary <- scenario_callwait_replication_summary %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished)) %>%
  mutate(scenario = "scenario")

# show top 10 results
head(scenario_callwait_summary, n = 10)

median_wait	mean_wait	var_wait	min_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countunfinished	countfinished	scenario
3.890618	3.289315	3.122619	0.001257	5	1.598594	5	0	1747.56	scenario

5.5 T-test

Use t.test() function to compare waiting times between baseline and scenario and see if there is a significant difference in the results.

5.5.1 Version taught in slidepack:

t.test(scenario_callwait_results$wait_time_call, baseline_callwait_results$wait_time_call)

## 
##  Welch Two Sample t-test
## 
## data:  scenario_callwait_results$wait_time_call and baseline_callwait_results$wait_time_call
## t = 16.071, df = 172339, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.1241180 0.1585979
## sample estimates:
## mean of x mean of y 
##  3.290622  3.149264

5.5.2 Paired T-test

This could actually be a paired t-test because they use the same common random number set and have the same number of replications therefore replication 1 of baseline and scenario are correlated (unless the scenario itself changes the andom number sampling). we should check the variances of the differences is less than the individual variances before using paired t test. To do it, we need to:

combine a dataset
compute the difference between the scenario and baseline mean waits for each pair of replications
conduct Shapiro-Wilk normality test for the differences and check p-value
If p-value is >0.05, we conduct paired t-test

part3_summary_byreplication <- bind_rows(
  baseline_callwait_replication_summary,
  scenario_callwait_replication_summary
)

d <- with(
  part3_summary_byreplication,
  mean_wait[scenario == "baseline"] - mean_wait[scenario == "scenario"]
)
shapiro.test(d) # => p-value = 0.6375

## 
##  Shapiro-Wilk normality test
## 
## data:  d
## W = 0.95276, p-value = 0.04435

res <- t.test(mean_wait ~ scenario, data = part3_summary_byreplication, paired = TRUE)
res

## 
##  Paired t-test
## 
## data:  mean_wait by scenario
## t = -8.5021, df = 49, p-value = 3.295e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1753380 -0.1082975
## sample estimates:
## mean of the differences 
##              -0.1418177

p-value is greater than the significance level 0.05 implying that the distribution of the differences (d) are not significantly different from normal distribution. therefore we do a paired T-Test

data:  mean_wait by scenario
t = -26.092, df = 49, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  -0.4228384 -0.3623634
sample estimates:
  mean of the differences
-0.3926009

The p-value of the test is < 2.2e-16, which is less than the significance level alpha = 0.05. We can then reject null hypothesis and conclude that the average wait in the baseline is significantly different from the average wait in the scenario

5.5.3 Boxplot of the call waits in the baseline vs scenario

Using the ggplot2 package and geom_boxplot(), create boxplots of the means of the replications of the call waiting time in the baseline compared to the scenario. Set theme and title if you wish. To compare boxplots, you can either

create 2 different plots and put them on one chart using gridExtra package

a <- ggplot(baseline_callwait_replication_summary) +
  geom_boxplot(aes(x = scenario, y = mean_wait),
               outlier.shape = 19,
               width = 0.3,
               fill = "white",
               colour = "orange") +
  labs(y = "mean call wait (mins)") +
  ylim(0, 4) +
  theme(axis.title.x = element_blank(),
        axis.text = element_text(size = 11))

b <- ggplot(scenario_callwait_replication_summary) +
  geom_boxplot(aes(x = scenario, y = mean_wait),
               outlier.shape = 19,
               width = 0.3,
               fill = "white",
               colour = "orange") +
  labs(y = "mean call wait (mins)") +
  ylim(0, 4) +
  theme(axis.title.x = element_blank(),
        axis.text = element_text(size = 11))

grid.arrange(a, b, ncol = 2, nrow = 1, top = "Call wait mean results of baseline(left) and scenario(right)")

create one plot from combined dataset. Using bind_rows(), create one dataset with both scenarios and plot them together using ggplot2 already have the averages of by replications in part3_summary_byreplication

part3_summary <- bind_rows(baseline_callwait_summary, scenario_callwait_summary)

# plot
g <- ggplot() +
  geom_boxplot(data = part3_summary_byreplication,
               aes(x = scenario, y = mean_wait),
               outlier.shape = 19,
               width = 0.3,
               fill = "white",
               colour = "orange") +
  geom_point(data = part3_summary,
             aes(x = scenario, y = mean_wait),
             shape = 19,
             size = 2,
             colour = "red") +
  labs(x = "scenario", y = "mean call wait (mins)", title = "Boxplot of mean call wait from 50 replications") +
  theme(axis.text = element_text(size = 11))

g

5.5.4 Part 3 Questions

Looking at the results of baseline of 50 replications, what can you say about the results?

What do you interpret from the results of the scenario in comparison to the baseline?

What other analysis might you want to do of the model results?

Additional analysis: example if you want to compare ECR to Phones waits might do it this way:

baseline_wait_results <- envs_baseline %>%
  get_mon_arrivals(ongoing = TRUE) %>%
  distinct() %>%
  full_join(baseline_attributes, by = c("name", "replication")) %>%
  filter(start_time >= 0) %>% # so use this to remove the not started e.g. start_time -1
  mutate(wait_time_call = case_when(
    reneged == 1 ~ 5,
    !finished ~ t_day(5) - start_time,
    TRUE      ~ pmax(end_time - start_time - activity_time, 0)
  ))

baseline_wait_replication_summary <- baseline_wait_results %>%
  group_by(replication, via_ecr) %>% #group by both of these to get results for phone and ECR
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(!finished),
            countfinished = sum(finished)) %>%
  mutate(scenario = "baseline")

## `summarise()` has grouped output by 'replication'. You can override using the `.groups` argument.

# show top 10 results
head(baseline_wait_replication_summary, n = 10)

replication	via_ecr	median_wait	mean_wait	var_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countunfinished	countfinished	scenario
1	0	3.825883	3.183036	3.479728	5.000	1.464796	5.00000	0	1726	baseline
1	1	56.571437	121.432726	60680.123013	1129.364	14.716072	99.92010	3	354	baseline
2	0	3.744117	3.221769	3.258414	5.000	1.555424	5.00000	0	1770	baseline
2	1	65.218754	149.913667	69329.798636	1080.604	15.965540	149.12634	1	291	baseline
3	0	3.775021	3.116034	3.682359	5.000	1.206570	5.00000	0	1772	baseline
3	1	32.526169	64.839190	22473.421606	1075.412	5.002446	64.15030	6	314	baseline
4	0	3.952666	3.235613	3.446565	5.000	1.535425	5.00000	0	1722	baseline
4	1	41.584256	83.401127	35084.985697	1099.012	12.369562	73.55387	1	303	baseline
5	0	3.810451	3.201863	3.433046	5.000	1.446567	5.00000	0	1736	baseline
5	1	40.110651	76.963237	27211.750428	1008.028	9.223911	84.78485	0	299	baseline

baseline_wait_summary <- baseline_wait_replication_summary %>%
  group_by(via_ecr) %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished))

baseline_wait_summary

via_ecr	median_wait	mean_wait	var_wait	min_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countunfinished	countfinished
0	3.765179	3.147497	3.580269	0	5.000	1.334405	5.00000	0.00	1722.80
1	42.015194	94.007905	41412.497885	0	1052.733	10.070932	85.24734	4.08	302.48

What other ways might you want to present data?

For example, plot results from all replications

plot(get_mon_arrivals(envs_baseline), metric = "waiting_time")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 1062 rows containing non-finite values (stat_smooth).

## Warning: Removed 1062 row(s) containing missing values (geom_path).

plot(get_mon_arrivals(envs_scenario), metric = "waiting_time")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 29002 rows containing non-finite values (stat_smooth).

## Warning: Removed 29002 row(s) containing missing values (geom_path).

Histogram

ggplot(part3_summary_byreplication, aes(mean_wait)) +
  geom_histogram(fill = "orange") +
  labs(title = "Distribution of average call wait in 50 replications") +
  facet_wrap(vars(scenario), ncol = 2)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

5.5.5 Additional questions if you want to explore model further

What other opportunities are there for improvement of the process?

What other scenarios might you want to consider with the practice? What would your performance indicators be? How would you incorporate them in the model? Have a go if you have time. warm up period, staff rota, percentage coming via ecr, include preemption and restart

What else would need to be considered before they could implement changes suggested by the model? costs, other tasks the staff do, feasibility of suggested rota in real life

Are we running the model for a sufficient number of runs? process converted to R from: Simulations: The Practice of Model Development and Use, Robinson, S., 2nd Ed. 2014, Palgrave Macmillan

Remember back to our t test results we saved in res gave a confidence interval based on student’s t distribution

res$conf.int[1]

## [1] -0.175338

res$conf.int[2]

## [1] -0.1082975

Create a blank data frame to hold the cumulative results

check_baseline_replications <- data.frame(matrix(ncol = 5, nrow = 0))
x <- c("replication", "cmean", "tL", "tU", "percentdeviation")
colnames(check_baseline_replications) <- x

#for loop to see how the percent deviation of the confidence limits from the mean vary by the number of replications
for (i in 2:50) {
  df <- baseline_callwait_replication_summary %>%
    dplyr::select(replication, mean_wait) %>%
    filter(replication <= i) %>%
    summarise(replication = i,
              cmean = mean(mean_wait),
              tL = t.test(mean_wait)$conf.int[1],
              tU = t.test(mean_wait)$conf.int[2],
              percentdeviation = (cmean - tL) / cmean * 100) %>%
    data.frame()
  check_baseline_replications <- bind_rows(check_baseline_replications, df)
}

# show top 10 results
head(check_baseline_replications, n = 10)

replication	cmean	tL	tU	percentdeviation
2	3.202403	2.956327	3.448479	7.684104
3	3.173613	3.040728	3.306499	4.187198
4	3.189113	3.103887	3.274340	2.672417
5	3.191663	3.133636	3.249691	1.818099
6	3.168606	3.094869	3.242343	2.327118
7	3.145290	3.062986	3.227594	2.616748
8	3.144801	3.075910	3.213691	2.190611
9	3.147488	3.087915	3.207060	1.892702
10	3.139632	3.084423	3.194840	1.758443
11	3.153205	3.095464	3.210947	1.831200

6 Part 3 if you use now() function

To include the set_attribute using now() we reference the correct environment for each replication

We can now analyse results

6.0.1 Count reneged

envs_baselineB %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  group_by(replication) %>%
  summarise(countreneged = n()) %>%
  ungroup() %>%
  summarise(countreneged_mean = mean(countreneged))

countreneged_mean
579.86

n=580

6.0.2 Count F2F

envs_baselineB %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  group_by(replication) %>%
  summarise(countF2F = n()) %>%
  ungroup() %>%
  summarise(countF2F_mean = mean(countF2F))

countF2F_mean
143.04

n=143

6.0.3 Call waits

baselineB_callwait_results <- envs_baselineB %>%
  get_mon_attributes() %>%
  pivot_wider(-c(time), names_from = key, values_from = value) %>%
  filter(via_ecr == 0, arrival >= 0) %>%
  mutate(wait_time_call = case_when(
    !is.na(reneged) ~ 5,
    is.na(call_taken) ~ (10 * 60) - arrival, # in queue/progress still, model run for 10 hours
    TRUE ~ call_taken - arrival # has been answered
  ))

baselineB_callwait_replication_summary <- baselineB_callwait_results %>%
  group_by(replication) %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(is.na(call_end) & is.na(reneged)),
            countfinished = sum(!is.na(call_end) | !is.na(reneged)))

# show top 10 results
head(baselineB_callwait_replication_summary, n = 10)

replication	median_wait	mean_wait	var_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countfinished
1	3.825883	3.183036	3.479728	5	1.464796	5	1726
2	3.744117	3.221769	3.258414	5	1.555424	5	1770
3	3.775021	3.116034	3.682359	5	1.206570	5	1772
4	3.952666	3.235613	3.446565	5	1.535425	5	1722
5	3.810451	3.201863	3.433046	5	1.446567	5	1736
6	3.550199	3.053321	3.688564	5	1.141632	5	1699
7	3.525437	3.005393	3.874786	5	1.005038	5	1652
8	3.808467	3.141373	3.616938	5	1.237069	5	1733
9	3.892609	3.168985	3.635974	5	1.225592	5	1750
10	3.623343	3.068928	3.714392	5	1.107349	5	1662

baselineB_callwait_summary <- baselineB_callwait_replication_summary %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished)) %>%
  mutate(scenario = "baseline") # to identify it later

baselineB_callwait_summary

median_wait	mean_wait	var_wait	min_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countunfinished	countfinished	scenario
3.765179	3.147497	3.580269	0	5	1.334405	5	0	1722.8	baseline

baselineB_callwait_summary$mean_wait

## [1] 3.147497

baseline_callwait_summary$mean_wait

## [1] 3.147497

6.1 Part 3 Scenarios with option B

Trajectory must be included inside the lapply for replication

We can now analyse results

6.1.1 Count reneged

envs_scenarioB %>%
  get_mon_attributes() %>%
  filter(key == "reneged") %>%
  group_by(replication) %>%
  summarise(countreneged = n()) %>%
  ungroup() %>%
  summarise(countreneged_mean = mean(countreneged))

countreneged_mean
603.16

n=752

6.1.2 Count F2F

envs_scenarioB %>%
  get_mon_attributes() %>%
  filter(key == "booked_appt", value == 2) %>%
  group_by(replication) %>%
  summarise(countF2F = n()) %>%
  ungroup() %>%
  summarise(countF2F_mean = mean(countF2F))

countF2F_mean
171.56

n=165

6.1.3 Call waits

scenarioB_callwait_results <- envs_scenarioB %>%
  get_mon_attributes() %>%
  pivot_wider(-c(time), names_from = key, values_from = value) %>%
  filter(via_ecr == 0, arrival >= 0) %>%
  mutate(wait_time_call = case_when(
    !is.na(reneged) ~ 5,
    is.na(call_taken) ~ (10 * 60) - arrival, # in queue/progress still, model run for 10 hours
    TRUE ~ call_taken - arrival # has been answered
  ))

scenarioB_callwait_replication_summary <- scenarioB_callwait_results %>%
  group_by(replication) %>%
  summarise(median_wait = median(wait_time_call),
            mean_wait = mean(wait_time_call),
            var_wait = var(wait_time_call),
            min_wait = min(wait_time_call),
            max_wait = max(wait_time_call),
            Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
            Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
            countunfinished = sum(is.na(call_end) & is.na(reneged)),
            countfinished = sum(!is.na(call_end) | !is.na(reneged)))

# show top 10 results
head(scenarioB_callwait_replication_summary, n = 10)

replication	median_wait	mean_wait	var_wait	min_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countfinished
1	3.749845	3.209906	3.181770	0.0000000	5	1.464999	5	1666
2	3.932084	3.337168	3.053690	0.0000000	5	1.716950	5	1770
3	3.890046	3.282838	3.134920	0.0000000	5	1.540661	5	1789
4	3.764143	3.236900	3.134823	0.0000000	5	1.548852	5	1712
5	3.901525	3.279950	3.099232	0.0000000	5	1.651958	5	1723
6	4.060991	3.371249	3.040204	0.0111140	5	1.751085	5	1764
7	4.017161	3.352447	3.025200	0.0009612	5	1.774429	5	1763
8	3.904756	3.298098	3.105040	0.0000000	5	1.632506	5	1755
9	4.038012	3.353227	3.106330	0.0000000	5	1.716001	5	1764
10	3.684992	3.213686	3.210820	0.0000000	5	1.500257	5	1708

scenarioB_callwait_summary <- scenarioB_callwait_replication_summary %>%
  summarise(median_wait = mean(median_wait),
            mean_wait = mean(mean_wait),
            var_wait = mean(var_wait),
            min_wait = mean(min_wait),
            max_wait = mean(max_wait),
            Qtile25th_wait = mean(Qtile25th_wait),
            Qtile75th_wait = mean(Qtile75th_wait),
            countunfinished = mean(countunfinished),
            countfinished = mean(countfinished)) %>%
  mutate(scenario = "scenario") # to identify it later

scenarioB_callwait_summary

median_wait	mean_wait	var_wait	min_wait	max_wait	Qtile25th_wait	Qtile75th_wait	countunfinished	countfinished	scenario
3.890618	3.289315	3.122619	0.001257	5	1.598594	5	0	1747.56	scenario

scenarioB_callwait_summary$mean_wait

## [1] 3.289315

scenario_callwait_summary$mean_wait

## [1] 3.289315

Workshop for Simmer in R training

Stacey Croft, Anastasiia Zharinova and Tom Jemmett

8th July 2021

1 Building the model - part 1

1.1 Libraries needed

1.2 Create distribution functions

1.3 Creating the trajectories for the branches

1.3.1 ECR branch

1.3.2 Phone branch

1.3.3 Booking option branches

1.4 Patient Trajectory

1.5 Add resources and generator to the environment

1.6 Run the model

1.7 Glimpse results to check run

1.8 Questions

1.8.1 Count reneged

1.8.2 Count F2F

1.8.3 Call waits

2 Scenario additional model features

3 Extra model features - part 2

3.1 New environment

3.2 Add time-dependent interarrival times

3.3 Phone branch with prioritisation

3.4 Patient Trajectory

3.5 Schedule for staffing

3.6 Add resources and generator to the environment

3.7 Run model

3.8 Questions

3.8.1 Count reneged

3.8.2 Count F2F

3.8.3 Call waits

3.9 What else would you do to improve this model?

3.10 What other metrics would you want to measure / the practice find useful?

3.11 Alternative way to calculate waiting time - option B

3.11.1 Part 2 Questions for Option B

4 Multiple runs and Scenarios - part 3

4.1 Replication of baseline

4.2 get baseline results and plots

4.2.1 Count reneged

4.2.2 Count F2F

4.2.3 Call waits

5 Scenario

5.1 create patient flow 3 for this scenario with 40% ECR

5.2 scenario interarrivals

5.3 lapply as before but using patient_flow3 and increasing frquency of arrivals

5.4 Scenario results

5.4.1 Count reneged

5.4.2 Count F2F

5.4.3 Call waits

5.5 T-test

5.5.1 Version taught in slidepack:

5.5.2 Paired T-test

5.5.3 Boxplot of the call waits in the baseline vs scenario

5.5.4 Part 3 Questions

5.5.5 Additional questions if you want to explore model further

6 Part 3 if you use now() function

6.0.1 Count reneged

6.0.2 Count F2F

6.0.3 Call waits

6.1 Part 3 Scenarios with option B

6.1.1 Count reneged

6.1.2 Count F2F

6.1.3 Call waits

8^th July 2021